# 1.ANOVA检验
setwd("C:\\Users\\lenovo\\Desktop\\R语言")    #设置文件路径
ADdata=read.csv("merdata.csv")      #读取ADdata数据
rownames(ADdata)=ADdata[,1]          #(以下两行)将第一列置为列名
ADdata=ADdata[-1]
ADdata=log2(ADdata)                  #取对数
ADdata[ADdata=="-Inf"]=NA            #将负无穷大值置于NA
group<-c()                     
for(i in 1:64)                    #(以下循环)读取列名并通过列名分组
{
  group[i]=substr(colnames(ADdata)[i],1,2)
}
ADdata=rbind(group,ADdata)       #合并(后续看来好像没啥必要)
group1=grep("as",ADdata[1,])      #(以下三行)确定三组分属那几列
group2=grep("ad",ADdata[1,])
group3=grep("ct",ADdata[1,])
j1=15                    #置值
k1=39
m1=10
p<-c()
for(i in 2:101)         #循环
{
  for(j in 1:15)        #3个小循环主要计算有数据点数
  {
    if(is.na(ADdata[i,group1[j]])) 
      j1=j1-1
  }
  for(k in 1:39)  
  {
    if(is.na(ADdata[i,group2[k]]))
      k1=k1-1
  }
  for(m in 1:10)
  {
    if(is.na(ADdata[i,group3[m]]))
      m1=m1-1
  }
  if(!(j1>2&k1>2&m1>2))     
    p[i-1]=NA       #单组非0数据点少于3个的蛋白，P值输出为NA
  else             
  {
    total=as.double(unlist(ADdata[i,]))    #total值
    group4=as.factor(ADdata[1,])         #group组
    res=oneway.test(total~group4)        #oneway.test检验    
    p[i-1]=res$p.value                  #得到P值
  }
  j1=15                     #(以下三行)循环一次后置初值
  k1=39
  m1=10
}
#最终P值向量即为所得
ADdata1=read.csv("merdata.csv")      #读取原始ADdata数据
p=rbind(as.vector(ADdata1[,1]),p)     #加入蛋白名一列
p=t(p)        #倒置



# 2.GO富集分析
setwd("C:\\Users\\lenovo\\Desktop\\R语言")    #设置文件路径
load("class5_volcano.RData")            #载入数据
data<-prostat[with(prostat,(P<0.05)),]       #筛选P<0.05的蛋白
if (!requireNamespace("BiocManager", quietly = TRUE))     #安装富集分析的包
  install.packages("BiocManager")
BiocManager::install("clusterProfiler")
BiocManager::install("org.Hs.eg.db")
library("clusterProfiler")
library("org.Hs.eg.db")
install.packages("enrichplot")
library("enrichplot")
install.packages("ggplot2")
library("ggplot2")

id=as.vector(data[,1])       #提出蛋白质ID那一列
entrezID=mget(id, org.Hs.egSYMBOL2EG, ifnotfound=NA)   #寻找基因号，未找到置NA
entrezID <- as.character(entrezID)
out=cbind(data,entrezIDs=entrezID)     #加入基因号那一列
out2<-c()   
for(i in 1:353)       #进行循环去除NA的蛋白
{
  if(!(out[i,4]=="NA"))
  {
    out2=rbind(out2,out[i,]) 
  }
}
gene=out2$entrezIDs       #获取最终基因名
kk <- enrichGO(gene=gene,        #GO富集分析
               
               OrgDb = org.Hs.eg.db,
               
               pvalueCutoff =0.05,
               
               qvalueCutoff = 0.05,
               
               ont="all",
               
               readable =T)
pdf(file="气泡图.pdf",width = 10,height = 8)      #保存为气泡图
dotplot(kk,showCategory = 10,split="ONTOLOGY") + facet_grid(ONTOLOGY~., scale='free')
dev.off()